#include "Scalar.h"


CScalar::CScalar(void)
{
	m_pScalar = NULL;
	m_w = m_h = m_d = 0;
}


CScalar::~CScalar(void)
{
	SafeDeleteArray(m_pScalar);
	m_w = m_h = m_d = 0;
}


int CScalar::OpenScalarFile(const char *pathname, FILE *fp)
{
	if (pathname == NULL)
		return -1;			// no input pathname

	fprintf(fp, "Reading Scalar file %s ... ", pathname);

	Nrrd *scalar = nrrdNew();
	if (nrrdLoad(scalar, pathname, NULL)) {
		fprintf(fp, "cannot load file %s:\n%s\n", pathname, biffGetDone(NRRD));
		return -2;			// cannot load file
	}

	m_w = (uint)(scalar->axis[0].size);
	m_h = (uint)(scalar->axis[1].size);
	m_d = (uint)(scalar->axis[2].size);

	m_pScalar = new float[m_w*m_h*m_d];
	memcpy(m_pScalar, scalar->data, sizeof(float)*m_w*m_h*m_d);
	nrrdNuke(scalar);

	fprintf(fp, "Success.\n");
	fprintf(fp, "Width: %d, Height: %d, Depth: %d\n", m_w, m_h, m_d);

	return 0;
}

int CScalar::SaveScalarFile(const char *pathname, FILE *fp)
{
	if (pathname == NULL || m_pScalar == NULL)
		return -1;

	Nrrd *dst = nrrdNew();
	size_t ss[3];
	ss[0] = m_w;
	ss[1] = m_h;
	ss[2] = m_d;

	nrrdAlloc_nva(dst, nrrdTypeFloat, 3, ss);
	float *data = (float *)(dst->data);
	uint size = m_w * m_h * m_d;
	memcpy(data, m_pScalar, sizeof(float)*size);

	nrrdAxisInfoSet_va(dst, nrrdAxisInfoSpacing, 1.0, 1.0, 1.0);
	nrrdAxisInfoSet_va(dst, nrrdAxisInfoKind, nrrdKindSpace, nrrdKindSpace, nrrdKindSpace);
	nrrdAxisInfoSet_va(dst, nrrdAxisInfoUnits, "mm", "mm", "mm");

	if (nrrdSave(pathname, dst, NULL)) {
		fprintf(fp, "cannot save file %s:\n%s\n", pathname, biffGetDone(NRRD));
		return -1;
	}

	nrrdNuke(dst);

	return 0;
}


float CScalar::GetScalar(const uint x, const uint y, const uint z) const
{
	if (x < 0 || x > (m_w-1) ||
		y < 0 || y > (m_h-1) ||
		z < 0 || z > (m_d-1)) {
		return 0.0f;
	}

	uint index = (z * m_h + y) * m_w + x;
	return m_pScalar[index];
}

float CScalar::GetScalar(const float x, const float y, const float z) const
{
	if (x < 0.0 || x > (float)(m_w-1) ||
		y < 0.0 || y > (float)(m_h-1) ||
		z < 0.0 || z > (float)(m_d-1)) {
		return 0.0f;
	}

	uint x1 = (uint)floor(x);
	uint y1 = (uint)floor(y);
	uint z1 = (uint)floor(z);

	uint x2 = (uint)ceil(x);
	uint y2 = (uint)ceil(y);
	uint z2 = (uint)ceil(z);

	float a = x - (float)x1;
	float b = y - (float)y1;
	float c = z - (float)z1;

	uint i1 = (z1 * m_h + y1) * m_w + x1;
	uint i2 = (z1 * m_h + y1) * m_w + x2;
	uint i3 = (z1 * m_h + y2) * m_w + x2;
	uint i4 = (z1 * m_h + y2) * m_w + x1;

	uint i5 = (z2 * m_h + y1) * m_w + x1;
	uint i6 = (z2 * m_h + y1) * m_w + x2;
	uint i7 = (z2 * m_h + y2) * m_w + x2;
	uint i8 = (z2 * m_h + y2) * m_w + x1;

	float v1 = m_pScalar[i1];
	float v2 = m_pScalar[i2];
	float v3 = m_pScalar[i3];
	float v4 = m_pScalar[i4];
	float v5 = m_pScalar[i5];
	float v6 = m_pScalar[i6];
	float v7 = m_pScalar[i7];
	float v8 = m_pScalar[i8];

	float scalar = v1 + a * (v2 - v1) + b * (v4 - v1) + c * (v5 - v1) + 
		a * b * (v1 - v2 + v3 - v4) + 
		a * c * (v1 - v2 + v6 - v5) + 
		b * c * (v1 - v4 + v8 - v5) + 
		a * b * c * (-v1 + v2 - v3 + v4 + v5- v6 + v7 - v8);

	return scalar;
}

float CScalar::GetScalar(const CPoint3F &p) const
{
	return GetScalar(p.m_x, p.m_y, p.m_z);
}

float CScalar::GetScalar(const CPoint3F *p) const
{
	return GetScalar(p->m_x, p->m_y, p->m_z);
}

int CScalar::CreateScalar(const uint w, const uint h, const uint d,
						  const float *pData)
{
	if (m_pScalar != NULL)
		return -1;

	uint size = w * h * d;
	if (size == 0)
		return -2;

	m_pScalar = new float[size];
	memset(m_pScalar, 0, sizeof(float)*size);
	m_w = w;
	m_h = h;
	m_d = d;

	if (pData != NULL) {
		memcpy(m_pScalar,pData, sizeof(float)*size);
	}

	return 0;
}

void CScalar::SetScalar(const uint x, const uint y, const uint z, const float s)
{
	if (m_pScalar == NULL)
		return;

	if (x < 0 || x > (m_w-1) ||
		y < 0 || y > (m_h-1) ||
		z < 0 || z > (m_d-1)) {
		return;
	}

	uint index = (z * m_h + y) * m_w + x;
	m_pScalar[index] = s;
}

void CScalar::SetScalar(const uint index, const float s)
{
	if (m_pScalar == NULL)
		return;

	m_pScalar[index] = s;
}

void CScalar::Normalize()
{
	uint size = m_w * m_h * m_d;
	float min_v = m_pScalar[0];
	float max_v = m_pScalar[0];
	for (uint i = 1; i < size; ++i) {
		if (m_pScalar[i] > max_v)
			max_v = m_pScalar[i];
		if (m_pScalar[i] < min_v)
			min_v = m_pScalar[i];
	}

	for (uint i = 0; i < size; ++i) {
		m_pScalar[i] = (m_pScalar[i] - min_v) / (max_v - min_v);
	}
}